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1 Abstract 

We report research experience in applying an Unsteady Reynolds-Averaged Navier-Stokes (URANS) solver for the prediction 
of time-dependent flows in the presence of an active flow control device. The configuration under consideration is a synthetic 
jet created by a single diaphragm piezoelectric actuator in quiescent air. Time-averaged and instantaneous data for this case 
were obtained at Langley Research Center, using multiple measurement techniques. Computational results for this case using 
one-equation Spalart-Allmaras and two-equation Menter’s turbulence models are presented here along with comparisons with 
the experimental data. The effect of grid refinement, preconditioning and time-step variation are also examined. 

2 Introduction 

Significant interest has been growing in the aerospace community in the field of flow control in recent years. An entire AIAA 
conference is now devoted every other year to this field. In particular, an international workshop, CFDVAL2004 [1], was held in 
March 2004 at NASA Langley Research Center to assess the state-of-the-art for measuring and computing aerodynamic flows 
in the presence of synthetic jets. Thomas et al. [2] have conducted an exhaustive and comprehensive survey identifying the 
feasibility of using active flow control to improve the performance of both external and internal flows. Suggested applications 
cover a wide range from smart materials and micro-electro-mechanical systems (MEMS) to synthetic zero net mass jets for 
enhancing control forces, reducing drag, increasing lift and enhancing mixing to name a few. It is also conjectured that active 
flow control would permit the use of thicker wing sections in non-conventional configurations, such as the Blended Wing Body 
(BWB) configuration without compromising the aerodynamic performance. 

Most of the research in the area of active flow control is of empirical nature, mainly due to the cost and lack of confidence 
in computational methods for such complex flows. However, without the availability of efficient and well-calibrated computa- 
tional tools, it will be a very difficult, expensive and slow process to determine the optimum layout and placement for active 
flow control devices in practical applications. With the continuous reduction of computer costs in recent years, more attention 
is being devoted to the simulation of such unsteady flows, and many researchers have recently started examining various flow 
control devices from a computational point of view (Refs. [3] - [8]). With few exceptions, most of the numerical studies are 
undertaken without an active interaction with experimental investigators. Comparisons with experimental data are sometimes 
done years after the experimental data have been acquired. Under such a scenario, one has to reconstruct some of the details 
about the experimental arrangement and boundary conditions without the benefit of concrete and consistent information. Based 
on our experience from previous validation exercises [9], we recognized the need for active collaboration of the computational 
and experimental research. Without a symbiotic relationship among such groups, major misunderstandings can develop when 

* Senior Member 

t Professor, Department of Mathematics, Associate Fellow 


1 



results from these disciplines indicate significant differences. We were very fortunate to have a frank and cooperative relation- 
ship with the researchers conducting the experiments as well as access to pertinent experimental data. 

Our primary objective for this work is to calibrate an existing computational scheme with experimental data for time- 
dependent flows encountered in active flow control environments. Special attention is devoted to establish appropriate and 
stable boundary conditions for such flows, especially in the absence of detailed experimental data required for closure. 

The configuration chosen for CFD validation is identified as Case 1 in the CFDVAL2004 workshop [1], and represents 
an isolated synthetic jet formed by a single diaphragm piezoelectric actuator exhausting into ambient quiescent air. Multiple 
measurement techniques including PIV, LDV and hot-wire probes were used to generate a large body of experimental data 
for this configuration. The details of the experimental setup and geometric configuration are described in Refs. [1, 10]. In 
this paper, we assess the effects of grid refinement, preconditioning and turbulence models on the flow field generated by this 
synthetic jet flow control device. We replace the actuator cavity with a simpler configuration. We demonstrate and calibrate our 
computational method for simulating synthetic jets by comparing the numerical results with the experimental data. 


3 Governing Equations 


A generalized form of the thin-layer Navier-Stokes equations is used to model the flow. The equation set is obtained from 
the complete Navier-Stokes equations by retaining the viscous diffusion terms normal to the solid surfaces in every coordinate 
direction. For a body-fitted coordinate system (£, rj, () fixed in time, these equations can be written in the conservative form as: 


„ ,3(U) , 3(F-F V ) , 3(G — G v ) , 3(H - H v ) 
Vd ^T + m + 3, + 3C = 


( 1 ) 


where U represents the conserved variable vector. The vectors F, G, H, and F v , G v , H v represent the convective and diffusive 
fluxes in the three transformed coordinate directions, respectively. In Eqn. (1), Vol represents the cell-volume or the Jacobian of 
the coordinate transformation. A multigrid-based, general purpose multi-block structured grid approach is used for the solution 
of the governing equations. In particular, the TLNS3D flow code is used in this study to solve Eqn. (1). Several references 
exist describing the discretization and algorithmic details of TLNS3D code. We include only a brief summary of the general 
features, and refer to the work of Vatsa and co-workers [11, 12] for further details regarding the TLNS3D code. 


4 Spatial Discretization 

The spatial terms in Eqn. (1) are discretized using a cell-centered finite volume scheme. The convection terms are discretized 
using second-order central differences with matrix artificial dissipation (second- and fourth- difference dissipation) added to 
suppress the odd-even decoupling and oscillations in the vicinity of shock waves and stagnation points [13, 14, 15]. The viscous 
terms are discretized with second-order accurate central difference formulas [11]. The zero-equation model of Baldwin-Lomax 
[16], one-equation model of Spalart-Allmaras [17] and Menter’s two-equation SST model [18] are available in TLNS3D code 
for simulating turbulent flows. For the present computations, the Spalart-Allmaras (SA) model and the Menter’s SST model are 
used for simulating turbulent flows. 


5 Temporal Discretization 

Regrouping the terms on right hand side into convective and diffusive terms, Eqn. (1) can be rewritten as: 

^ = -C(IJ) + D P (U) + D a (U) (2) 

at 

where C(U), D p (U), and D a ( U) are the convection, physical diffusion, and artificial diffusion terms, respectively. The 
cell-volume or the Jacobian of the coordinate transformation is included in these terms. 

The time-derivative term can be approximated to any desired order of accuracy by a Taylor series 

^ = ^1 k> u " +1 + oiU" + a 2 \J n - 1 + a 3 U"- 2 + ...] (3) 
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The superscript n represents the last time level at which the solution is known, and n+1 refers to the next time level to 
which the solution will be advanced. Similarly, n-1 refers to the solution at one time level before the current solution. Eqn. (3) 
represents a generalized backward difference scheme (BDF) in time, where the order of accuracy is determined by the choice of 
coefficients ao, a 1 , a 2 ... etc. For example, ao = 1.5, a\ = —2 and 02 = .5, results in a second order accurate scheme (BDF2) 
in time, which is the primary scheme chosen for this work due to its robustness and stability properties [19]. Regrouping the 
time-dependent terms and the original steady-state operator leads to the equation: 


^ TT „ +1 ^(U"’"- 1 -) 

At At 


5(U”) 


(4) 


where E(U n,n 1 •) depends only on the solution vector at time levels n and earlier. S represents the steady state operator or 
the right hand side of Eqn. (2). By adding a pseudo-time term, we can rewrite the above equation as: 


m , ^0 + i ffCU"-"- 1 -) 

dr At At 


S(U n ) 


(5) 


6 Solution Algorithm 

The algorithm for solving unsteady flow relies on the steady-state algorithm in the TFNS3D code [11, 12], The basic algorithm 
consists of a five-stage Runge-Kutta time-stepping scheme for advancing the solution in pseudo-time, until the solution con- 
verges to a steady state. Efficiency of this algorithm is enhanced through the use of local time-stepping, residual smoothing and 
multigrid techniques developed for solving steady-state equations. Because the Mach number in much of the domain is very 
low we consider the use of preconditioning methods [26, 27]. 

In order to solve the time-dependent Navier-Stokes equations (Eqn. 5), we add another iteration loop in physical time 
outside the pseudo-time iteration loop in TLNS3D. For fixed values of S(U”), E( U"’” -1 ’"), we iterate on U n+1 using the 
standard multigrid procedure of TLNS3D developed for steady-state, until the desired level of convergence is achieved. This 
strategy, originally proposed by Jameson [20] for Euler equations and adapted for the TLNS3D viscous flow solver by Melson 
et. al [19], is popularly known as the dual time-stepping scheme for solving unsteady flows. The process is repeated until 
the desired number of time-steps have been completed. The details of the TLNS3D flow code for solving unsteady flows are 
available in Refs. [19, 25, 28]. 


7 Boundary Conditions 

The boundary conditions required for solving the Navier-Stokes equations, such as the no-slip, no injection, fixed wall tem- 
perature or adiabatic wall, far-field and extrapolation conditions are well understood and readily available in most flow codes 
including the TLNS3D code. The most accurate procedure to simulate the moving diaphragm would require moving grid ca- 
pability. For simplicity, we chose to simulate this type of boundary condition by a periodic velocity transpiration condition. 
The frequency of the transpiration velocity at the diaphragm surface in the numerical simulation corresponds to the frequency 
of the oscillating diaphragm. The peak velocity at the diaphragm surface was obtained from numerical iteration to match the 
experimentally measured peak velocity of the synthetic jet emanating from the slot exit. The pressure at the moving diaphragm 
is also required for closure. However, in the absence of unsteady pressure data from the experiment, we imposed a zero pressure 
gradient at the diaphragm boundary. We also tested the pressure gradient boundary condition obtained from one-dimensional 
normal momentum equation [21], which had very little impact on the solutions. Due to the simplicity and robustness, we 
selected the zero pressure gradient boundary condition at the diaphragm surface. 

8 Synthetic Jet Test Case: Background 

The test configuration examined in this paper corresponds to a single diaphragm piezoelectric actuator operating in quiescent air. 
The oscillatory motion of the diaphragm produces a synthetic jet that exhausts into surrounding air. This configuration, shown 
in Fig. 1, consists of a 1.2 mm wide rectangular slot connected to a cavity with a piezoelectric diaphragm, and corresponds 
to case 1 of the CFDVAL2004 workshop on flow control devices [1], The cavity and diaphragm geometry of this actuator 
are highly three-dimensional in the interior. However, the actual slot through which the fluid emerges is a high aspect ratio 
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rectangular slot and can be modeled as a two-dimensional configuration. Partial view of the 2-D grid provided to the workshop 
participants, and used by the present authors is shown in Figs. 2(a) and 2(b). The computational results contributed by the 
workshop participants are available in refs. [1, 22]. A consensus developed during the workshop that simulating the flow 
field inside the actuator cavity with an oscillating piezoelectric diaphragm from first principles was beyond the capability of 
the existing CFD codes. Most of the workshop participants, including the present authors modeled the internal cavity of the 
actuator as a two-dimensional configuration, and simulated the diaphragm motion via a transpiration condition imposed at 
the diaphragm surface. Some of the workshop participants further simplified the cavity modeling by imposing a transpiration 
condition at the bottom part of the slot’s neck or even directly at the slot exit. After examining these results, we concluded 
that as long as the unsteady velocity signal at the slot exit replicates experimental conditions, details of the cavity modeling 
have an insignificant effect on the development of the synthetic jet emanating from the slot. Another conclusion derived from 
the CFDVAL2004 workshop was that no particular methodology or turbulence model emerged superior for simulating this test 
case [1, 22], 

One of the major difficulties identified during the CFDVAL2004 workshop was the large disparity in experimental data 
obtained from PIV, hot-wire probes and LDV measurement techniques. Such a variation in experimental data made it difficult 
to validate the numerical methods. Part of the difficulty in acquiring a consistent set of experimental data arose from the fact that 
the performance of the piezoelectric diaphragm depends on ambient conditions. Therefore, its performance degrades over time 
which means that for a given input voltage, the actuator produces smaller jet velocities as it ages. Because these experiments 
were conducted over a period of several months, inconsistencies creeped in the data. 

We include here sample results obtained with the TLNS3D code to encapsulate the status of CFD simulations for this 
configuration at the conclusion of the CFDVAL2004 workshop. For these computations, we used the 2-D grid of approximately 
65,000 nodes from the CFDVAL2004 workshop website as the baseline grid. As depicted in Figs. 2(a), 2(b), this grid includes 
the internal cavity and the diaphragm geometry. The transpiration condition is imposed at the diaphragm surface, which is 
represented by the vertical line on the left side at the bottom of Fig. 2(b). We present the time-history of the (vertical) v- 
velocity at x = 0, y — 0.12 mm, and the time average of the v-velocity along the jet centerline in figures 3 and 4, respectively. 
The computations were performed with several grid densities, two physical time steps and two turbulence models. The coarse 
grid (eg) was obtained by eliminating every other point from the baseline grid, and a fine grid (fg) was obtained by adding 50% 
points in the normal direction. We see from these figures that varying all these factors did not produce any major changes in 
the computational results in region near the slot exit. However, further away from the slot exit, coarse grid (eg) results indicate 
that the coarser grid does not provide adequate resolution for this problem. A finer grid (fg) and a reduction in time step (low 
dt) had an insignificant effect on the computational results. Because of the scatter between the experimental PIV and hotwire 
data, it is difficult to make any definitive conclusion regarding the accuracy of turbulence models. 

9 Results 

Yao et. al [10] have recently revisited the synthetic jet test case and acquired experimental data for this configuration with a new 
diaphragm. The detailed field data was obtained with the PIV technique. In addition, pointwise data along the jet centerline 
was obtained with hotwire and LDV techniques. The performance of the actuator was monitored regularly. The results from 
this study show good consistency among multiple measurement techniques [10]. 

We simulated the new experimental test case with a simplified cavity geometry, shown schematically in Fig. 5. The tran- 
spiration condition is imposed at the bottom of the slot’s neck to simulate the velocity generated by the oscillating diaphragm. 
Similar boundary condition treatment produced satisfactory results at the CFDVAL2004 workshop, and has also been studied 
in detail by Yamaleev and Carpenter [23]. They demonstrated that for actuators with deep cavities, specifying the transpiration 
condition at a distance of at least 4-5 slot widths away from the slot exit produces only a small loss in numerical accuracy. 
A top-hat velocity profile, with a dominant frequency of 450 Hz., replicating the experimental conditions was imposed at the 
bottom boundary. The precise form of the velocity signal was obtained by curve fitting the measured velocities at the slot exit 
(x=0, y=0.3 mm) with a Fast Fourier Transform to reflect the proper mode shape and to ensure zero net mass transfer. The 
amplitude of this velocity was determined numerically to match the peak velocity from the experiment at the slot exit. The 
free stream Mach number in the exterior quiescent region is specified as M x = .001 for simulating incompressible flow in the 
compressible flow code to avoid numerical difficulties at Mach zero. 

The computational grid for this case was obtained from the baseline grid described in the previous section by eliminating 
the portions of grid below the neck of the slot. The resulting grid consists of over 60,000 nodes, and should provide adequate 
resolution based on the grid refinement study reported in reference [1], and as seen from the results in figures 3 and 4. Similarly, 
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72 time-steps/period corresponding to 5° phase angle between the time steps should provide adequate temporal resolution for 
this problem. The one-equation Spalart-Allmaras (SA) [17] turbulence model was used to simulate the effect of turbulence in 
these computations. Based on the peak jet velocity and slot width, the Reynolds number is approximately 3000, and therefore 
falls in the regime where the jet is expected to be turbulent. Therefore, the flow was assumed to be fully turbulent in the present 
simulations. 

The time-history of the vertical velocity for a complete period from the computational results with SA model is compared 
with the experimental data in Fig. 6 at x=0 and y=0.3 mm. This is the closest point to the slot exit where the PIV data is 
available. In addition to the PIV data, LDV measurements are also available at this location and are shown here. The LDV 
data was scaled down by a factor of 0.9 as suggested by Yao et. al [10] to match the diaphragm displacement for the two 
sets of measurements, and is in good agreement with the PIV data except for a flatter region near the 150° phase. The overall 
agreement between the computational and experimental results is quite good at this location, and therefore serves as an accurate 
boundary data for the jet emerging from the slot. 

Next, we compare the time-averaged v-velocities along the jet centerline in Fig. 7 with the PIV and LDV data. The 
experimental data from two different techniques (PIV and LDV) is found to be in good agreement with each other, giving 
credibility to the accuracy and consistency of the measurements for this case. The overall agreement of the baseline TLNS3D 
results with the experimental data is also quite good. The low-speed preconditioning results included in this figure are essentially 
identical to the baseline results. Because low-speed preconditioning ([26]-[28]) primarily reduces the artificial viscosity for 
unsteady flows, it is inferred that the artificial viscosity is already low in these simulations. The time-averaged v-velocity 
profiles at y=l and 4 mm are shown in Fig. 8 and Fig. 9, respectively. Except for a smaller velocity peak at the centerline at 
y=l mm, the computational results are in very good agreement with the experimental data. 

The contour plots of the time-averaged v-velocities based on PIV measurements covering a distance of 8 mm from slot exit 
are shown in Fig. 10(a). Note that the high resolution PIV data was limited up to a distance of 8 mm from slot exit because 
of cost and time constraints. The contour plots from our baseline computations are shown in Fig. 10(b) for comparison. The 
computational results appear to accurately capture all of the prominent features seen in the PIV data including the width and 
spreading rate of the synthetic jet. 

We now examine the phase-averaged velocities at selected locations in space and time, starting with v-velocities at y=2 and 
4 mm along the jet center-line. The PIV and LDV data along with baseline TLNS3D solutions are shown at these locations in 
Figs. 1 1(a) and 11(b). The computational results are in fairly good agreement with the two sets of experimental data, especially 
in the suction phase. The agreement with the experimental data further away from the slot exit is slightly worse during the 
peak expulsion cycle. In particular, the CFD results predict a delayed phase shift for the peak expulsion, reflective of smaller 
convective speed for outward movement of the vortex compared to the experimental data. 

We gain a broader perspective of the flow-field by examining the contour plots of the velocities at the phase angles represen- 
tative of the expulsion (phase=75°) and suction (phase=255°) cycles. The velocity contours for streamwise (u-vel) and vertical 
velocities (v-vel) obtained from the PIV data and TLNS3D computational results are shown in Figs. 12 - 15 and Figs. 16 - 19. 
These figures were generated using identical contour levels for both the experimental and CFD data for providing a quantitative 
comparisons. Note that solid lines represent positive values, while dashed line represent negative values for the velocities. This 
sign convention is helpful in identifying the flow direction and the position of the vortex center. It is clear from these figures 
that the computational results capture most of the pertinent features observed experimentally and are in very good agreement 
for the suction phase. During the expulsion phase, the computed vortex center is located closer to the slot exit compared to 
experimental data, although the peak velocity at the vortex center is in good agreement with the PIV data. Yao et al. [10] have 
observed increasing three-dimensional effects for this case as one moves away from the slot exit, mainly due to ring vortices 
formed from the slot edges. We conjecture that these ring vortices induce forces that accelerate the convection of synthetic jet 
in the far field. Wheteher this indeed is the primary reason for differences in convective speed of the vortex, can only be verified 
by 3-D simulations with sufficient resolution to capture the vortical structure emanating from the edges of the slot. 


10 Concluding Remarks 

Detailed comparisons have been presented for time-averaged and phase-averaged velocities between experimental data and 
CFD results. The effect of truncation errors were found to be small based on a grid refinement, preconditioning and physical 
time-step refinement studies. The differences between the solutions obtained from the one -equation turbulence model of Spalart 
and two-equation model of Menter were found to be small in the near field for the synthetic jet. The modeling of the internal 
flow in the cavity of the actuator turned out be an extremely difficult problem. Fortunately, the development of the synthetic 
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jet in the quiescent medium is driven primarily by the velocity field at the slot exit, and detailed modeling of the cavity is not 
warranted. Based on comparisons of the computational results with the original experimental data, it was decided to repeat 
the experiments to obtain more consistent data. The computational results in the reduced domain with a modified forcing 
function are found to be in much better agreement with the new experimental data in the near field. However, the agreement 
with the experimental data deteriorates in regions further away from the slot exit. Based on the available experimental data, it 
appears that the flow becomes three-dimensional after 5-6 slot widths away from the exit. Future work should focus on 3-D 
computations for this configuration to resolve such issues. 
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Figure 1 : Schematic of Piezoelectric Actuator 
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(a) Global view (b) Detailed view 

Figure 2: Computational grid for Piezoelectric Actuator 
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Figure 3: Time-history of v-velocity at x=0, y=0.12 Figure 4: Time-averaged v-velocity along jet center- 
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Figure 5: Simplified model of actuator 
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Figure 6: Time-history of v-vel near slot exit. Synjet II 



Figure 8: Average v-vel at y=l mm. Synjet II 



Figure 7: Average v-vel along centerline. Synjet II 



Figure 9: Average v-vel at y=4 mm, Synjet II 
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